Stochastic wave function approach to generalized master equations 



Heinz-Peter Breuer, Bernd Kappler and Francesco Petruccione 
Albert-Ludwigs-Universitat, Fakultat fur Physik, 
Hermann-Herder Strafle 3, D-79104 Freiburg ira Breisgau, Federal Republic of Germany 

(February 1, 2008) 

A generalization of the stochastic wave function method is presented which allows the unravelling 
of arbitrary linear quantum master equations which are not necessarily in Lindblad form and, 
moreover, the explicit treatment of memory effects by employing the time-convolutionless projection 
operator technique. The crucial point of this construction is the description of the open system in 
a doubled Hilbert space, which has already been successfully used for the computation of multitime 
correlation functions. 
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I. INTRODUCTION 

Usually, the state of an open quantum system is de- 
scribed by a reduced density matrix p(t) which is a pos- 
itive operator on the Hilbert space Hs of the system. 
On the other hand, within the stochastic wave function 
method the state of the open system is described by an 
ensemble of pure, normalized states ip(t) the covariance 
matrix of which equals the reduced density matrix [Q-^) , 



P (t)= / DipDr\Tp)^\Pbi>,t}. 



(1) 



In Eq. (Q) the integral extends over the Hilbert space of 
the system, DipDip* denotes the Hilbert space volume 
element, and P[ijj,t] is the time-dependent probability 
density of finding the state of the system in the volume 
element DiftDip* near ip M. This formulation has essen- 
tially two advantages compared to the conventional de- 
scription: First, this approach allows the investigation of 
the dynamics of an individual quantum system which is 
continuously observed by some measurement device , 
whereas the reduced density matrix can only describe the 
state of an ensembles of quantum systems. Second, from 
a computational point of view, the numerical integration 
of the quantum master equation can become rather ex- 
pensive for large systems, since the reduced density ma- 
trix has ./V 2 degrees of freedom, where N is the dimension 
of the system's Hilbert space. In contrast, a stochastic 
wave function has only N components, which can signif- 
icantly reduce the computational expense J8j. Moreover, 
algorithms which are based on stochastic simulations can 
easily be implemented on parallel computers. 

The dynamics of the stochastic wave function is gov- 
erned by a stochastic evolution equation, and the con- 
struction of this evolution equation within the Born- 
Markov approximation is well understood. However, 
in some situations non-Markovian effects can signifi- 
cantly alter the reduced system dynamics. In this ar- 
ticle we will present a scheme which allows a system- 
atic incorporation of memory effects into the stochas- 
tic wave function method. To this end, we make use 
of an expansion scheme which is known from the the- 



ory of non-equilibrium statistical mechanics - the time- 
convolutionless projection operator technique [p|Jl0|l. 

This article is organized as follows. In Sec. || we discuss 
the unravelling of quantum master equations by stochas- 
tic wave functions. This concept is well known for quan- 
tum master equations which are in Lindblad form fl I 
and we brie fly summarize the major results in Sec. |lA. 
In Sec. [I B we generalize this concept to the treatment 
of arbitrary linear quantum master equations. Using this 
result, we present in Sec. Ill a general framework which 
allows an explicit treatment of memory effects within 
the stochastic wave function method. This concept is 
then illustrated by means of an exactly solvable model in 
Sec. IV - the damped Jaynes-Cummings model. 



II. STOCHASTIC SIMULATION OF QUANTUM 
MASTER EQUATIONS 

A. Quantum master equations in Lindblad form 

In the Markovian regime, the time evolution of the re- 
duced density matrix p(t) is governed by the quantum 
master equation in Lindblad form 



at'® 



(2) 



+ !>*(*) \-\^\U P (t) - \pit)L\L l + L t p{t)A , 
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where Hq is the Hamiltoninan of the system, the time- 
dependent coefficients Si(t) describe an energy shift in- 
duced by the coupling to the environment, namely the 
Lamb and Stark shifts, and the positive rates ji(t) model 
the dissipative coupling to the i— th decay channel. This 
evolution equation is either obtained by a phenomeno- 
logical ansatz or through a derivation which is based on 
a microscopic model of the system-reservoir interaction. 

Using similar techniques, one can also obtain a Marko- 
vian time evolution equation for the stochastic state vec- 
tor ip. There are several phenomenological approaches 
which simply construct a stochastic evolution equation 



1 



in such a way, that the equation of motion of the co- 
variance matrix is the quantum master equation 
This procedure is often called unravelling of the quan- 
tum master equation ffl . Other approaches are based on 
a continuous observation of the system under considera- 
tion by some measurement device, for example a photon 
detector, and employ the basic measurement postulates 
for the description of the dynamics of an individual quan- 
tum system HQ . Finally, similar to the derivation of the 
quantum master equation, it is also possible to obtain 
the stochastic time evolution directly from an underly- 
ing microscopic model by an explicit derivation of the 
differential Chapman Kolmogorov equation for the prob- 
ability density P[ip,t] [§. 

A particular example of such a stochastic evolution 
equation which arises in the above approaches is the 
stochastic differential equation 



d^{t) = -iG(i/j,t)dt + J2 



il>{t) dNi{t) 7 



(3) 



where the diV,*(t) are the differentials of indepen- 
dent Poisson process Ni(t) with mean (dNi(t)) — 
j^WL^fdt and 

Gty, t) = H(t)i/> + S ^) L \ L ^ 

i 



This particular equation of motion describes the time 
evolution of a piecewise deterministic process. The dif- 
ferential of the Poisson process dNi(t) can either take 
the value or 1. If diVj(t) = 0, then the system evolves 
continuously according to the nonlinear Schrodinger-type 
equation 



(5) 



whereas, if dNi(t) = 1 for some i, then the system un- 
dergoes an instantaneous, discontinuous transition of the 
form 



U^{t) 

w^mw 



(6) 



Note that the generator G(tp, t) of the continuous time 
evolution is non-Hermitian and hence the propagator of 
ip(t) is non-unitary. However, due to the nonlincarity of 
the generator, the norm of ip{t) is preserved in time. 

Using the standard Ito calculus for the differentials 
dNi(t) of a Poisson process, i. e., dNi(t)dNj(t) = 
6ijdNi(t), it is easy to check, that the equation of motion 
of the covariance matrix of ip(t) equals the usual Marko- 
vian quantum master equation for the reduced density 



matrix in Lindblad form. Thus, both descriptions yield 
the same equations of motion for the expectation values 
of system observables. Finally, we want to remark that 
it is also possible to extend the stochastic wave functions 
method to the calculation of arbitrary matrix elements 
of system operators and hence to the determination of 
multitime correlation functions p|,|l2]]. 

B. General quantum master equations 

In this section we present a generalization of the 
stochastic wave function method to quantum master 
equations which are not in Lindblad form (in Sec. Ill 
we will also encounter this type of evolution equations). 
To be more specific, we consider an equation of motion 
for the reduced density matrix of the form 



^p(t)=A(t)p(t)+p(t)B^t) 



J2c t (t)p(t)D\(t), (7) 



with some arbitrary time-dependent linear operators 
A(t), B(t), Ci(t), and Di(t). This form represents the 
most general linear equation of motion for p(t), which 
is local in time, i. e., an equation of motion where p(t) 
only depends on p(t). In order to find an unravelling 
of this equation of motion we follow a strategy, which 
has already been successfully applied to the calculation 
of multitime correlation functions : We describe the 
state of the open system by a pair of stochastic state 
vectors 9 — (</>,-0) T which is an element of the doubled 
Hilbert space H = H © H, in such a way, that 



P (t) = / D0D9*\(/>)(iP\P[e,t\, 



(8) 



where the integral extends over the doubled Hilbert space 
H, and P[9, t] is the probability density of finding the sys- 
tem in the "state" 9 at time t. Furthermore, we define 
the operators F(t) and Ji(t) as 



_ / A(t) 
B{t) 



Ji{t) 



Ci(t) 
Di{t) 



(9) 



An unravelling of the quantum master equation (]?]) by a 
stochastic wave function 9(t) can be obtained using the 
stochastic differential equation 



d6{t) = -iG(0,t)dt 



(10) 



|0(*)ll 



\\Mt)9(t) 



-Ji(t)0(t) -9{t) dN.it) 



where dNi (t) is the differential of a Poisson process with 
mean 



(dNiit)) = 



\Mt)9(t)\ 



-dt, 



(11) 
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and 



G(e,t)=i[F(t) + lYl 



\\Mt)Q(tW 
\m)\\ 2 



9{t). (12) 



Again, the stochastic differential equation ( fTo|) describes 
a piecewise deterministic jump process, where G(9, t) is 
the generator of the continuous time evolution and the 
operators Ji (t) lead to discontinuous instantaneous tran- 
sitions. In order to show that the stochastic differen- 
tial equation (|l^) leads to the correct equation of motion 
for p(t) one can rewrite Eq. ( [Tof ) as a system of coupled 
stochastic differential equations for </> and tp, and com- 
pute the mean of the differential 



d {[</>) = + \4>)(dip\ + |#)<#| 



(13) 



using the Ito calculus, which justifies our ansatz. It is 
important to note that this unravelling also contains the 
unravelling presented in Sec. II A . If the equation of mo- 
tion for the reduced density matrix is in Lindblad form, 
and (f>(0) = ip(0) then it is easy to see that <f)(t) and 
ip(t) are identical for all t and the stochastic differential 
equation ([lO]) reduces to Eq. (^) with the choice 



A(t) = B(t) = -iH s - \ hk(t) + iS k {t)} h\h k 



and 



Ci(t) = Di(t) 



(14) 



(15) 



III. NON-MARKOVIAN STOCHASTIC WAVE 
FUNCTION METHOD 

In this section we present a general scheme which 
allows a systematic, perturbative treatment of mem- 
ory effects within the stochastic wave function method. 
This approach is based on the time-convolutionless pro- 
jection operator technique p|jlQ[|, which is related to 
the Nakajima-Zwanzig projection operator technique 

mi- 

As a microscopic model, we consider a system which 
is coupled to an environment. The Hamiltonian of the 
total system is given by 



H = H Q + aHj, 



(16) 



where H describes the free evolution of the system and 
the reservoir and H\ their interaction. The parameter a 
denotes a dimensionless expansion parameter. The state 
of the total system is described by the interaction picture 
density matrix W(t) which is a solution of the Liouville- 
von Neumann equation 



d_ 



W(t) = -ia[Hi(t),W(t)] = aL(t)W(t), (17) 



where the interaction Hamiltonian in the interaction 
picture is defined as Hi(t) = exp(iH t)Hiexp(—iHot). 
Since we are interested in the dynamics of the reduced 
system, we define a projector V as 



VW{t) = Tr R {W{t)} <g> pr = p{t) ® pr, 



(18) 



where pr is a stationary state of the environment, and 
a projector Q — 1 — V . A quasi-closed equation of mo- 
tion for VW(t) can be obtained by using the Nakajima- 
Zwanzig projection operator technique fl3| [l5[], namely 



dsK{t,s)PW(s) 



where the memory kernel K(t 1 s) is defined as 

K(t,s) = a 2 PL(t)g(t,s)L(s), 
with the propagator 



Q (t, s) = T^ exp 



ds'QL(s') 



(19) 



(20) 



(21) 



and T^ indicates the chronological time ordering. In ob- 
taining Eq. (|l9|) we have assumed that TrR{i/j 2fc+1 pR} = 
0, and that the system and the reservoir are uncorrelated 
initially, i. e., WJO) = p(0) ® Or. To eliminate the time- 
convolution PJlCjl in Eq. ([l9]) we replace W(s) by the 
expression 

W(s) =G(t,s)(T + Q)W(t), (22) 
where the backward propagator G(t, s) is defined as 



G{t, s) = T^ exp 



a / ds'L(s') 



(23) 



and T_^ indicates the anti-chronological time ordering. 
This leads to the time-convolutionless equation of mo- 
tion 



§-VW(t) = K(t)VW(t), 

where the generator K(t) is defined as 
K(t) = aVL(t)[l-Z(t)}- 1 V 

and 

E(t)=af dsg{t,s)L{s)VG(t,s) 
Jq 



(24) 



(25) 



(26) 



If the operator (1 — can be expanded in a geo- 

metric series, which is possible if the coupling between 
the system and the reservoir is not too strong, then we 
can rewrite the generator K(t) as 
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K{t) = aY J VL{t){Y.{t)) n V (27) 
and obtain a perturbative expansion in the form 

oc 

K{t) = Y. a2nK ^)- ( 28 ) 

n=0 

Note that all terms containing odd orders of the coupling 
constants vanish, since by definition of V and Lit) we 
have VLit\) ■ ■ ■ Lit 2k +i)P = 0. The explicit expressions 
for the second and fourth order contribution arc 



K 2 {t) = j dbxPL{t)L{ti)V, 
Jo 



and 



Ki{t) = dti dti \ dt 3 



(29) 



(30) 



JO JO 

-PL(t)L(t 2 )PL{h)L(t a )P - VL{t)L{t z yPL{tx)L{t2)fP 



The higher order contributions can be obtained in a sys- 
tematic way by a slight modification of van Kampen's 
cumulant expansion jl6| (see also ||). 

The time-convolutionless quantum master equation 
( p4|) allows us to use the stochastic wave function method 
for the description of the dynamics of the open system. 
To this end, we note that the equation of motion for the 
reduced density matrix pit) which results from either us- 
ing Eq. ( p4| ) directly or any perturbative approximation 
of this equation, is linear in p(t) and local in time. Hence, 
we can write it in the form of Eq. (|t]) where, of course, 
the operators A{t), Bit), and Di(t) depend on the 

interaction Hamiltonian Hi, and apply the unravelling of 



this equation of motion described in Sec. [IB. This leads 



to a stochastic wave function description of the open sys- 
tem which can be formulated, at least in principal, to 
any desired order in the coupling, and which is hence 
generally applicable to any open quantum system. 



IV. EXAMPLE 

As a specific example for the general concept presented 
in Sees, [il] and III we consider the spontaneous decay of 
a two-level system coupled to the electromagnetic field 
which is initially in the vacuum state within the rotat- 
ing wave approximation. The Hamiltonian of the total 
system is given by 



Ho = u) S (T + o +y^ y Ukbjb k , (31) 

k 

Hi = <j+ <g) B + <t~ ® B f with B = ^g k b k , (32) 



where cos and 0J k denote the eigenfrequencies of the sys- 
tem and reservoir, respectively, and the g k are real cou- 
pling constants. As usual, denote the pseudospin op- 
erators, and the b k are the annihilation operators for the 
field mode k. Inserting the above definitions into the ex- 



pressions for K 2 {t) and K^t), Eqs. ( j29| ) and (30), we ob- 
tain an equation of motion for the reduced density matrix 
pit), the time-convolutionless quantum master equation 



p(t) = - l -SW(t)[a + a-,pit)] 



d_ 

dt 

+ 7 (4) W 



(33) 



-cr + (7 pit) pit)a + a +a p(t)a~* 



which is in Lindblad form with time-dependent coeffi- 
cients. The coefficients to fourth order are given by 



«fti*(t-ti) + - 
o 1 



SW(t) 

x hit - t 2 Mh - t 3 ) + *(t - t a )*(*i - *a) 



dti / dt 2 / dt% 
ii Jo Jo 



+ * (i - t 8 )*(ti - i a ) + *(* - t 3 )*(ti - t 2 ) (34) 

and 

7 (4) (i) = / <fti$(t-*i) + ^ / dh 



1 



*1 /"*2 

dt2 / dt-< 



o * Jo Jo Jo 

*(* - *a)*(*i - h) - *(* - *a)$(*i - h) 

*(* - t 3 )*(ti - t 2 ) - *(t - t 3 )$(t! - t 2 ) 



(35) 



The real functions $(i) and ^(t) are related to the reser- 
voir correlation functions through 



$(t) + i¥(t) = 2Tr R {S(t)fit pR } e -st 
2 / ^J(w)e l(ws - W) *, 



(36) 



where -B(£) = cxpiiHot)B cxp(— iH$t), and we have per- 
formed the continuum limit. The function J(w) is the 
spectral density times the strength of the coupling at 
the frequency uo. This function determines the statistical 
properties of the system dynamics. 

For this particular model an exact equation of motion 
for the reduced density matrix can be obtained in the 
following way: First we define the three states [jl7| 

V^o = |0) s ® |0) R 

1>i = |i)s®|o) R 

Tpk = |0) s ® |fc>R (37) 

where |0}s and |l)g indicate the ground and excited state 
of the system, respectively, the state |0)r denotes the vac- 
uum state of the reservoir, and |/c)r = &I|0)r denotes the 
state with one photon in mode fc. If an initial pure state 



4 



4>(0) can be expanded in terms of these states, then the 
state at time t has the form (l7| 



<f>(t) = CqV'O + Cx(t)lpi + ^ Ck(t)lpk, 



(38) 



with some probability amplitudes cq, ci(t), and Ck{t) and, 
hence, the reduced density matrix is given by 



Pit) 



|ci(t)| : 



Cl (t)c* 



(39) 



Differentiating p(t) with respect to time leads to a quan- 
tum master equation in the Lindblad form 



(40) 



j t p{t) = - l -S{t)[*+o-,p{t)] 

+7(*) \-\° + °~p{t) - \p{t)(T + a-+a-p{t)a+ 



where the time-dependent coefficients S(t) and j(t) are 
given by 

ff ( t ) = -29/^1, 7(t ) = -2R (41) 



Ci(i) 



d(t) 



It is important to note that the time-convolutionless ex- 
pansion of the equation of motion (^3|) reproduces the 
structure of the exact equation of motion ( ftcj ) to all 
orders in the coupling. On the other hand, a pertur- 
bative expansion of the exact Nakajima-Zwanzig equa- 
tion to fourth order (l^] also contains terms of the form 

(T + 0~ p(s)<7 + (J~ . 

In order to find a stochastic unravelling of the time- 
convolutionless quantum master equation, we have to 
distinguish two cases: if the function J (to) leads to a 
rate "f^ 4 \t) which is positive for all t, then we can use 



the "usual" stochastic unravelling described in Sec. [I A 



Otherwise, if the rate also takes negative values, we have 
to use the more general algorithm which we presented in 
Sec. [IB. We will illustrate both cases by means of the 



damped Jaynes-Cummings model, which describes the 
coupling of a two-level system to a single cavity mode, 
which in turn is coupled to an environment. For this 
model, the function J(uj) is given by 



i 



7oA 2 



A 2 ' 



2ir (u>o — lu) 2 

where ujq is the center frequency of the cavity. 
Eq. ( p6| ) we obtain 

$(t) = 7o Ae _At cos(A£), 
*(£) = 7 Ae _At sin(Ai), 



(42) 
Using 

(43) 
(44) 



where the detuning A is defined as A = los — ojq. Obvi- 
ously, for A = 0, the function is purely exponential, 
and ty(t) vanishes. Hence, j^' is positive for all t. On 
the other hand, for A ^ the rate 7^(i) oscillates and 
can take negative values if A is sufficiently large. We will 
discuss both cases separately. 



A. Damped Jaynes-Cummings model on resonance 



For the resonant damped Jaynes-Cummings model, 
the rate j( 4 \t) can be calculated using Eq. (pa). It is 
given by 

7 (4) (*) = 70 {l - e- At + ^ [sinh(Ai) - At] e - At } , (45) 

which we have illustrated in Fig. |l| for A = 570 together 
with the rate ^ 2 \t) and the exact decay rate 



7(*) 



27oAsinh(di/2) 



dcosh(cft/2) + Asinh(dt/2) ' 



(46) 



where d = ^J\ 2 — 270A. 



The exact rate is obtained by 
into Eq. @. 



inserting the amplitude c\ (t) (see 17 

Since the rate 7*- 4 ^ (t) is positi ve fo r all t we may use 
the unravelling presented in Sec. II A. The dynamics of 



the stochastic wave function is governed by the stochastic 
differential equation (^), where for this model the gener- 
ator G{tj),t) of the deterministic motion is given by 



and the instantaneous transitions lead to jumps of the 
form 



\w-mw 



\0h 



(48) 



i. e., the state of the system is projected onto the ground 
state. In Fig. ^ we illustrate the deviation of pn(t) 
from the Markovian population exp(— 70 4) for an ini- 
tially excited system. Obviously, the perturbative expan- 
sion of the time-convolutionless generator K (t) converges 
rapidly and leads to an excellent agreement of the ex- 
act solution and the solution of the time-convolutionless 
quantum master equation to fourth order. In Fig. |^ we 
also show the solution of the stochastic simulation of 
Eq. (H) for 10 5 realizations, which is in very good agree- 
ment with the solution of the time-convolutionless quan- 
tum master equation (|33|). 



B. Damped Jaynes-Cummings model with detuning 



Making use of Eq. ( 35 ) we can calculate the decay rate 



7^ 4 '(i) for the damped Jaynes-Cummings model with de- 
tuning, which yields 



5 



7 «(t) 



7oA 2 



A 2 + A 2 



[1 



-A/ 



(cos(Ai) 



■ sin 



(At))] 



2(A 2 + A 2 ) 3 
2 

2 



{[l-3(f) 2 ]( 



e" At cos 



(2Ai)) 



1 - (f) 4 |Aicos(Ai) 



1+0$) Aisin(At) 



#)' 



e" At sin 



in(2At)). 



(49) 



In Fig. ||we have depicted the rate ^Ht) together with 
the exact decay rate. The parameters are chosen such 
that performing the usual Born-Markov approximation 
leads to the constant decay rate 7m = 1- For short times, 
the decay rate shows damped oscillations and converges 
in the long time limit to a time-independent decay rate 
which is close to 7m- 

Note however, that in this case the decay rate can 
also take negative values, and the corresponding quan- 
tum master equation has to be unraveled using the pro- 
cedure described in Sec. [IB. Thus, the dynamics of the 
stochastic wave function 6{t) = ((j>(t), ip(t)) T , being an 
element of the doubled Hilbert space TL — H (B H is gov- 



erned by the stochastic differential equation (10), where 
the operators F and J are given by 



V. SUMMARY 

In this article we have presented a generalization of the 
stochastic wave function method to arbitrary linear quan- 
tum master equations, which allows an explicit treatment 
of memory effects in a systematic way. This is done by 
employing the time-convolutionless projection operator 
technique, which yields a perturbative expansion of the 
equation of motion of the reduced density matrix. The 
latter is then unraveled by a stochastic wave function 
in the doubled Hilbert space. By means of the damped 
Jaynes-Cummings model, which is an exactly solvable 
model, we have illustrated the general theory and tested 
the performance of this method. 
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a + a~ 



and 



J 



_(l [i) {t)v 




(50) 



(51) 



Here, the generator of the deterministic motion is given 
by 



G(6,t) = i[F+- 



i \\jo(t)\\' 



2 \\9(t)\ 



9{t). 



(52) 



and the jumps induce instantaneous transitions of the 
form 



9(t) 



\\J0(t)\\ {) \ I0)s 



(53) 



Hence, if the rate "f^(t) is positive, the jump leads 
to a positive contribution to the ground state popula- 
tion poo{t), whereas a negative rate leads to a negative 
contribution to poo- The results of a stochastic simu- 
lation of Eq. ( px| ) with 10 5 realizations is displayed in 
Fig. [| together with the analytical solution of the time- 
convolutionless quantum master equation to fourth or- 
der (|3^) and the exact solution, which are in very good 
agreement. This clearly demonstrates the usefulness of 
the simulation algorithm presented here. 
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FIG. 1. Decay rate of the excited state in the damped fig. 3. Decay rate of the excited state in the 

Jaynes-Cummings model using the time-convolutionless mas- damped Jaynes-Cummings model with detuning using the 

ter equation to second (TCL 2) and fourth (TCL 4) order, time-convolutionless master equation to fourth order (TCL 

compared with the exact decay rate. The reservoir correla- 4) compared with the exact decay rate. The parameters are: 

tion time is r R = O.^ 1 . 7o = 65, A = 19.5, and A = 8A. 
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FIG. 2. Deviation of the population pn(t) from the Marko- 
vian population e~ 70 ' for the damped Jaynes-Cummings 
model for the inverse reservoir correlation time A = 570: 
stochastic simulation (symbols) and exact solution of the 
time-convolutionless master equation to second (TCL 2) and 
fourth (TCL 4) order, compared with the exact solution. 
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FIG. 4. Decay of the population pn(t) for the damped 
Jaynes-Cummings model with detuning using the time-convo- 
lutionless master equation to fourth order (TCL 4) compared 
with the exact decay rate. The parameters are: 70 = 65, 
A = 19.5, and A = 8A 
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